Black-hole binary simulations: the mass ratio 10:1 
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We present the first numerical simulations of an initially non-spinning black-hole binary with 
mass ratio as large as 10 : 1 in full general relativity. The binary completes approximately 3 orbits 
prior to merger and radiates (0.415 ± 0.017) % of the total energy and (12.48 ± 0.62) % of the 
initial angular momentum in the form of gravitational waves. The single black hole resulting from 
the merger acquires a kick of (66.7 ± 3.3) km/s relative to the original center of mass frame. The 
resulting gravitational waveforms are used to validate existing formulas for the recoil, final spin and 
radiated energy over a wider range of the symmetric mass ratio parameter n — M-lM?./{M\ + M2) 2 
than previously possible. The contributions of i > 2 multipoles are found to visibly influence the 
gravitational wave signal obtained at fixed inclination angles. 

PACS numbers: 04.25.D-, 04.25.dg, 04.30.Db 



I. INTRODUCTION 

Following the breakthroughs of 2005 [H, IS 9 , the nu- 
merical relativity community has constructed several in- 
dependent codes @,HH&la&[l3,[Il[Il and proceeded 
at a breathtaking speed in generating new insights into 
the dynamics of inspiralling and coalescing black-hole bi- 
naries (BBH) and the resulting gravitational wave pat- 
terns |13| . For example, numerical relativity has provided 
vital information reg arding the kick or recoil in astro- 
physical mergers [M 03 03 03 El 03 M, US Ipjp, 
simulations of spin flip and precession phenomena [la, [24| 
and the interpretation of g ravitational waveforms to be 
observed in the future [Hf US US ES ES H3- See also 
[3ll [32l [33l . [341 . [H[ for binary simulations involving neu- 
tron stars. This progress has come timely, as earthbound 
gravitational wave detectors LIGO, VIRGO, GEO600 
and TAMA US 13, US HI are now collecting data at or 
near the design sensitivity. The combined use of theoret- 
ical predictions and sophisticated data analysis methods 
will be crucial in achieving the first direct detection of 
gravitational waves and thus opening a new window to 
the universe. 

Purpose of this work is to extend the range of mass ra- 
tios probed by numerical relativity to q = M1/M2 = 10 
corresponding to 77 = q/(l + q) 2 = 0.0826. This mass 
ratio is of particular interest for a variety of reasons. 
First, studies of the supermassive black hole formation 
history starting with light seed black holes predict a sig- 
nificant fraction of mergers in the range 1 < q < 10 
and, depending on accretion details, the possibility of a 
peak near q = 5. ..10 Similarly, detailed statistical 

analysis of the mass distribution of supermassive galac- 
tic black holes predicts that most mergers will occur in 
the range 3 < q < 30 [4l|. Finally, accurate numeri- 
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cal simulations with q = 10 will facilitate unprecedented 
comparisons with approximative calculations based on 
post-Newtonian expansions (cf. [IS 5S 0, SH| as well 
as 46] for a review) and extreme mass ratio predictions 
based onperturbation theory and self force calculations 
(see [U SS US US US US [Hi and references therein). 
By considering, for example, an expansion in the mass 
ratio parameter, a naive estimate of the error in pertur- 
bative calculations would be of the order of 1/q 2 or 1 % 
relative to background quantities, or 10 % for quantities 
of first order in 1/q. While such comparisons are beyond 
the scope of this paper, we will lay the foundation for 
future work by giving a detailed convergence analysis of 
the numerical results including estimates for the uncer- 
tainties. 

We will also use our results to probe recently pub- 
lished formulas for calculating the final spin and recoil 
resulting from the coalescence of two black holes with 
given initial physical parameters. Spin measurements of 
black holes via direct astrophysical observations have so 
far pro vided information about several individual holes 
[541. l55l. US US US US HI but appear as yet to be in- 
sufficient for constructing reliable spin-distribution func- 
tions. The community has therefore pursued the alterna- 
tive path of using theoretical predictions in the context 
of t he g rowth history and accretion processes of the holes 
[iiL HE HI, HI, [HE H3 • It is important for the model- 
ing of the individual binary mergers in such simulations 
to have available mappings between the initial param- 
eters of the binary and the final spin and recoil of the 
post-merger remnant. A better understanding of the dis- 
tribution function of the black-hole recoil also generates 
a great deal of interest in its own right. In particular, 
there remain open questions as to how generically recoil 
velocities of thousands of km/s result from astrophysi- 
cally realistic binary mergers. Such large recoil would 
not only result in intergalactic populations of black holes 
but also affect the central structure of galaxies and put 
severe constraints on possible scenarios of the black hole 
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formation history [H, H, [70|, [zH, [zl, [zl Izl- For recent 
discussions on direct observational signatures of recoiling 
black holes, see also [zl, [H, [z3 . 

Existing formulas predicting the kick and final spin 
as functions of the initial parameters of the binary are 
based on simulations using mass ratios 1 < q < 4 or 
0.25 > rj > 0.16 and, in the case of Baker et al. [73], using 
also q = 6 or r/ = 0.122 for the final spin. Our results 
clarify the validity of these formulas in the case of initially 
non-spinning holes. We further analyze the multipolar 
structure of the resulting gravitational waveforms and 
illustrate the significance of the higher order multipolcs 
in the gravitational wave signal. 

The paper is organized as follows. We describe in 
Sec. |TT] the numerical simulations together with a cali- 
bration of their uncertainties. In Sec. IIII1 we analyze the 
gravitational waveforms generated in the inspiral with 
regard to the total amount of energy, linear and angular 
momentum radiated and test the validity of existing fit- 
ting formulas for final spin and recoil. We further demon- 
strate that a significant fraction of the energy is radiated 
in higher order multipoles which implies that gravita- 
tional wave signals differ significantly from the typical 
quadrupole shape presented in most of the numerical rel- 
ativity literature. We conclude in Sec. [W] and discuss 
possible future studies based on the data presented in 
this work. 



II. NUMERICAL SIMULATIONS 

We have performed our simulations using the BAM code 
as described in [!, [7§| . In order to obtain sufficient accu- 
racy for this demanding type of simulations, we use sixth 
order discretization for the spatial derivatives [8(| and 
fourth order accurate integration in time. Initial data 
are provided by the spectral solver described in Ref . [8l[ . 
Gravitational waves are calculated in the form of the 
Newman-Penrose scalar ^4 according to the procedure 
described in Sec. Ill in [i| (but see for an alternative 
method to extract ^4). 

The model we are focusing on in this work represents 
a black-hole binary with initial parameters as given in 
Table HI We follow the convention of [83| and normalize 
initial parameters relative to the total black hole mass 
M = Mi + Mi and dimensional diagnostic quantities 
by their total initial values; the Arnowitt-Deser-Misncr 
(ADM) mass Mabm and the total initial angular mo- 
mentum J; n i. The initial tangential linear momentum P 
has been calculated from Eq. (65) of [4[ which gives a 
3rd order post-Newtonian estimate for the momentum of 
a quasi-circular configuration. We calculate the number 
of orbits completed by this configuration from the wave- 
form as described in Sec. Ill C of [HI and obtain about 
3 orbits and 6 wave cycles. 

The mass ratio q = 10 does not require any changes in 
the construction of initial data. However, let us point out 
one important issue that affects the puncture evolution 



method. We use the "000" gauge advection choice, that is 
we evolve the shift according to do/3* = jB z and doB 1 = 

doT 1 — r] s B l with do — dt — fi % di- In second order form 
the shift condition is 

dip = -a,!* - Vs doP\ (i) 

from which it is immediate that the physical dimension 
of the shift damping parameter r/ s is 

bh] = l/M. (2) 

(We have added the label "s" to distinguish r/ s from the 
mass ratio rj used elsewhere.) The parameter r/ s was 
introduced to control the dynamics of the shift vector 
[84 , |85j; in particular it influences the degree of slice 
stretching that develops near the black holes during dy- 
namic gauge evolution, but r/ s also affects the drift of 
coordinates near the outer boundary. Only certain val- 
ues of r] s lead to stable evolutions. This also applies to 
evolutions in the moving puncture framework, see e.g. Q 
for a discussion of some of the T] s dependence found in 
our evolutions, and [86] for related discussions. 

The issue that arises for q = 10 is that rj s is chosen to 
be a global constant, say 77,, = 2.0/M, but the effect of 
■q s on the slice stretching near the black holes is different 
if the black hole masses are different. Assuming that 
rj s = 2.0/M is a good choice near the black holes for equal 
mass, Mi = M2 ~ M/2, then increasing Mi by a factor 
of 5 means that rj s has to be replaced by r/ s /5 to obtain 
the same amount of slice stretching. Similarly, if Mi is 
reduced, then 77^ should be correspondingly enlarged. 

For these reasons the standard choice of r\ s = 2.0/M 
made in [3] did not work for q = 10 evolutions; the ef- 
fective rj s near the black holes did not lead to stable and 
accurate evolutions. A tell-tale sign for 77,,. being too large 
is if the orbits drift outwards rather than spiral inwards, 
which is accompanied by a loss of convergence with time. 
Numerical experiments revealed that rj s can be chosen 
in a certain interval without loss of convergence, which 
warrants a detailed study in the future. For the present 
purpose it is sufficient to note that choosing 

r) s = 1.375/Af (3) 

works near both black holes even for q = 10. We plan to 
investigate possible benefits of a non-constant, position 
dependent r/ s (as already suggested in [Ii],[85j]) in a future 
publication. 

A computational issue arising for increasing mass ra- 
tios is that such simulations typically require more com- 
puter time than simulations for comparable masses. The 
resolution requirement and the time step size is deter- 
mined by the smaller mass, while the physical time re- 
quired for one orbit increases with the total mass. We 
illustrate this by giving rough estimates for the number 
of computational time steps required for simulating the 
last 3 orbits prior to merger. In the equal-mass case 
q = 1, reasonable accuracy can be obtained by using a 
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TABLE I: Initial parameters and main results of the 10:1 mass ratio simulation studied in this work. 7711,2 and Mi,2 are the 
bare masses and black hole masses respectively. The Bowen-York parameters for the tangential linear momentum P and the 
coordinate separation D are normalized with respect to the total black hole mass M = Mi ± M2. Radiated energy and angular 
momentum are normalized with respect to their total values for the spacetime. Finally, we give the dimensionless spin and kick 
parameter of the merged hole. 



777 i 7772 Mi M2 MadM 


P/M 


D/M 


Brad 

-Ma dm 


Jrad 
Jtot 


Jfin fkick 


2.4831 0.2303 2.5 0.25 2.7381 


0.0415 


7.0 


(0.415 ±0.017)% 


(12.48 ±0.62)% 


0.259 ± 0.003 66.7 ± 3.3 km/s 



resolution h — Mj 48 on the finest refinement level corre- 
sponding to about 25, 000 time steps for the last 3 orbits. 
For the present simulation with q = 10, on the other 
hand, we obtain the same number of orbits after about 
250, 000 time steps in the medium resolution case labeled 
N = 68 below. It is imperative, therefore, to maintain 
high numerical accuracy over a longer evolution time as 
we increase q. A similar increase in computational de- 
mands was noted by Lousto & Zlochower who use 
an additional refinement level for spinning binaries with 
<7 = 4 (see their Sec. III). 

For these reasons, it is not permissible to infer infor- 
mation on the numerical uncertainties of our simulations 
by merely extrapolating convergence studies of equal- 
or mildly unequal-mass binaries as obtained for exam- 
ple in [j, [25| . Instead, we need to study convergence 
of the present scenario using three resolutions. In the 
notation of Sec. VI A of 4], our grid setup is given by 
Xf7a=i-375p x iV : 6 x 2N : 6], where the number of grid 
points is N = 60, 68 and 76 for the low, medium and 
high resolution runs, respectively. In Fig. [1] we show the 
resulting convergence analysis obtained for the I = 2, 
m = 2 mode of the Newman-Penrose scalar and the ra- 
diated energy E and linear momentum P extracted at 
r cx = 36.5 JVf adm- For this purpose, we have decom- 
posed the wave signal according to 



tively of J; n i. These values are well modeled by a function 
flo + ai/Vex which gives us a relative error of about 3.5 % 
for the value extracted at the largest radius 36.5 Madm- 
Similar results are obtained for the other radiated quan- 
tities. We investigate the phase error of the 22 mode due 
to the extraction radius in the same way. We consider the 
phase as a function of retarded time u = t — r cx and fit for 
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t/M A , 



1P22 = Ae 



i(f> 



(4) 



and studied amplitude A and phase <fr separately. We 
observe 6 th order convergence for both quantities and 
the total radiated energy. For the linear momentum 

P = Px + Py > tne convergence is closer to 5 th order. 
We estimate the error due to discretization by comparing 
the high resolution solution with the Richardson extrapo- 
lated result assuming 4 th order convergence. We use this 
more conservative 4 th order extrapolation to account for 
4 th order accurate ingredients in the BAM code. The un- 
certainties thus obtained are about 5 % for the phase, 
1 % for the amplitude and radiated energy and 3 % for 
the recoil and radiated angular momentum. We empha- 
size that no alignment in phase or time of the wave signal 
has been applied for this analysis. 

We similarly determine the error arising from extract- 
ing the waves at finite radius. The extraction radii avail- 
able for this analysis are r ox = 18.3, 27.4 and 36.5 Madm- 
For the radiated angular momentum, for example, we ob- 
tain at these radii J rac i = 12.06, 12.31 and 12.46 % respec- 
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FIG. 1: Convergence plots for the amplitude A and phase <j> 
of the £ — 2, 777 = 2 mode of the Newman Penrose scalar 
(upper panels) as well as the radiated energy E and linear 
momentum P (lower panels). The scaling factors for fifth and 
sixth order convergence are Qs = 2.04 and Q@ = 2.30. 
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each value of u the function <p(u, r cx ) = 4>q(u)+<P\{u) / 'r ex - 
The difference between the extrapolated phase and the 
result obtained at r cx = 36.5 Madm is shown in Fig. [2] 
The accumulated phase error after coalescence is about 
A<p = 1.5 rad corresponding to a relative error of about 
2 %. The magnitude of these errors obtained at compar- 
atively small extraction radii seems to suggest that for 
q = 10 certain near zone effects are reduced compared to 
q = 1, although this is to be examined more closely in 
future work. 

We combine the uncertainty estimates due to dis- 
cretization and extraction radius assuming standard er- 
ror propagation and adding the squares of the individual 
errors. We thus obtain error estimates of about 4 % for 
the wave amplitude and radiated energy, 5 % for the ra- 
diated momenta and 5.5 % for the phase. Note that most 
of the phase error builds up during the late stages of the 
inspiral and merger, so that phase uncertainties for per- 
forming a comparison with post-Newtonian results will 
be smaller than the estimate given here. We summarize 
the results for the radiated energy, momenta and final 
spin m Table U with error bars based on the above anal- 
ysis. 



III. GRAVITATIONAL WAVE EMISSION 

In this section, we study in more detail the gravita- 
tional wave emission of the binary. We also discuss the 
resulting values for radiated energy, final spin and recoil 
in the context of general formulas suggested in the liter- 
ature. 



n 1 1 1 r 

- ^( r e x = 36 - 5 M AD U ) 




FIG. 2: Difference in phase of the t = 2, m = 2 mode 
extrapolated to infinite extraction radius and extracted at 
r cx = 36.5 Madm- 



A. Comparison with phenomenological formulas 

All theoretical modeling of the mass and spin evolu- 
tion of black holes in the context of their merger history 
requires a mapping between initial and final black hole 
parameters for each individual merger. The generation of 
such mappings has recently become an industrious area 
of research, largely because of the breakthroughs in nu- 
merical relativity which make it possible now to simulate 
black-hole binaries accurately through inspiral, merger 
and ringdown. Results are currently available for a rel- 
atively small subset of the parameter space only, and 
have resulted in various efforts to "extrapolate" to wider 
ranges of the inp ut p arameter space using (semi-)analytic 
methods [H, M, H HE M, M, HI • The fitting formulas 
thus generated, have been relatively well tested in the 
regime of binaries with nearly equal mass or in the test 
particle limit, but in the regime in between, correspond- 
ing to a symmetric mass ratio in the range?? = 0.05. ..0.12, 
accurate data have as yet not been available. Here, we 
fill this gap and test existing predictions for the case of 
a non-spinning binary with rj = 0.0826 or q = 10. 

In Fig. [3l we show the results for the energy radiated 
in the form of gravitational waves during the last three 
orbits of the inspiral, the merger and the ringdown as a 
function of 77. To be precise, we start the integration at a 
phase A(f> = 43.3 rad in the I = 2, m = 2 mode prior to 
the maximum amplitude in that mode. Similar results 
for mass ratios 1 < q < 4 were presented in [l4| and 
approximated with a polynomial fit in Eq. (3.13) in [25| . 
Our simulations with q = 10 are in excellent agreement 
with the amount of gravitational wave energy expected 
from the polynomial fit. We observe similar agreement 
with the extrapolated prediction of the peak of the energy 




FIG. 3: Predictions for the energy radiated during the last 
three orbits (starting at phase 43.3 radians before the maxi- 
mum in the I = 2, m = 2 mode), merger and ringdown are 
compared with numerical data from [l4| as well as new results 
obtained for the mass ratio q = 10 : 1. 
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flux given in Baker et al. [78|]. Using their relation (6), 
we obtain £ , 22,max = 3.26 x 10 -5 to be compared with our 
numerical result £^2,max = 3.18 X 10~ 5 . The difference 
is comfortably within either study's error estimates. 

In Fig. 2J we consider fitting formulas for the final spin 
of the merged hole taken from Berti et al. [H| , Rezzolla et 
al. Hi] , Tichy & Marronetti [H[ and Buonanno et al. (95j 
(Eq. (C6) in Baker et al. [73) as well as the EOB predic- 
tions by Damour & Nagar [96|]. The figure demonstrates 
that all formulas agree very well with both, the results of 
our previous study in [Tij in the range 77 = 0.16. ..0.25 and 
the new value at rj — 0.0826. On the other hand, we are 
not able to discriminate between the different formulas 
based on our results of non-spinning binaries. The figure 
also contains the predicted final spin of the model pro- 
posed by Buonanno, Kidder & Lehner (BKL) [94]. The 
values have been obtained by numerically solving their 
Eqs. (2-5). Given that they do not fit existing numerical 
data but model the final spin, their agreement with the 
numerical results is remarkable. Even for nearly equal 
masses, the deviations are relatively small, of the order 
of 5 %. 

Fig. E] shows a similar analysis for the gravitational re- 
coil using analytic predictions by Gonzalez et al. [Til ]. 
Baker et al. [93] and Schnittman & Buonanno [98| . We 
thereby additionally test Eq. (2) of Lousto & Zlochower 
23 1 which by construction reduces to the prediction of 
14j in the limit of initially non-spinning holes. We fur- 
ther emphasize that Schnittman & Buonanno model the 
recoil using the effective one body (EOB) method in- 
stead of merely fitting available numerical data (cf. BKL 
for the final spin above). It is natural in this case to 




Eqs.(2-5), BKL [94] 
Eq.(3.17b) Berti et al. [25] 
Tab.1, Damour, Nagar [96] 
Tab.1, Damour, Nagar [96] 
Eq.(5) Rezzolla et al. [89] 
Eq.(7) Tichy, Marronetti [91] 
Eq.(C6) Baker et al. [78] 

Gonzalez et al. [14] 

q=10 



FIG. 4: Same as Fig. [3] but for the final spin. The curves re- 
sulting from E q. (7) of Tichy & Marronetti and from Eq. (C6) 
of Baker et al. 78, 95] would be indistinguishable in this plot 
and have been represented as the single dash-dotted line. The 
same applies to the curves resulting from Berti et al [25[ and 
the upper predictions of Damour & Nagar's (9f| which as- 
sumes the Kepler law (see their Sec. II). 



expect, larger deviations from the numerically obtained 
values. In analogy to the comparison of the final spin, 
we observe generally good agreement between predicted 
values and the numerical results including q = 10. This 
agreement is encouraging as we believe it highly unlikely 
that there exist local extrema in either curve in the range 
rj = 0.05... 0.15 and therefore provides strong support 
for all formulas in the limit of vanishing initial spin of 
the holes. 



B. Multipolar structure 

Gravitational waves are commonly decomposed into 
multipoles using spherical harmonics of spin weight —2. 
Specifically, the complex Newman-Penrose scalar ^4 ex- 
tracted at constant radius r ex is written as a sum (see 
e. g. EH ) 



ipe m (t)- 2 Y im (d,(f)). 



(5) 



For illustration, we show in Fig. [5] the real part of the 
multipole coefficients for £ — m — 2,..., 5. Most of 
the earlier studies of black-hole binaries focused on non- 
spinning, equal-mass systems where > 98% is radiated 
in the quadrupole terms I = 2. For more general classes 
of binaries, however, a significant fraction of the gravi- 
tational wave energy is radiated in higher order modes, 
e.g. (25)]. Simultaneously, the gravitational wave signal 
predicted for a given orientation of the binary's orbital 
plane will contain substantial contributions from higher 
order multipoles and may thus exhibit a pattern much 
more complex than visible in the quadrupole amplitude 
^2±2(t). While this complex multipolar structure places 
higher demands on the modeling of the gravitational 
wave sources, it provides us with a large amount of infor- 
mation in the effort to estimate parameters from gravi- 
tational wave observations. 




FIG. 5: Same as Fig. [3] but for the recoil. 
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A comprehensive study of exploiting such information 
of higher order multipoles in the context of gravitational 
wave data analysis is beyond the scope of this paper. In 
this section we therefore restrict ourselves to a discussion 
of the multipolar distribution of the gravitational wave 
energy and merely illustrate the significance of including 
higher order modes in the waveform for fixed inclination 
angle of the orbital plane relative to an observer. 

The amount of energy contained in different multipoles 
for the inspiral of binaries with mass ratio q — 1...4 has 
been given in Table IV of [25[. 

We graphically display the energy contained in the 
multipoles in Fig. [7] Here, the upper panel displays the 
energy contained in the multipoles corresponding to a 
particular value of I, whereas the symbols in the lower 
panel show the fractional energy contained in all modes 
up to i < An ax . We discard the first 50 Madm of the 
waveform which is dominated by spurious radiation in- 
herent to the initial data. This corresponds to the ra- 
diated energy starting at phase 43.3 rad prior to the 



maximum in the I = 2, m = 2 mode. The resulting 
energy is well approximated by quadratic polynomials. 
Specifically, we obtain the following fits for the numer- 
ical data in the range 1 < q < 10 corresponding to 
0.25 > rj > 0.0826 
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- 49.97? - 
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Ee<5/E ra d 
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The trend for higher order multipoles to carry a larger 
fraction of the total radiated energy £" ra d is clearly main- 
tained for q = 10. Close inspection of Eg = ^ reveals a 
local maximum near 7? = 0.04. We believe this to be 
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FIG. 6: Multipole coefficients r ex MAnMipem for I = m 
2, ...,5 extracted at 36.5 Madm- 



X 


/ =2 

max 


+ 


1 =3 

max 


O 


1 =4 

max 


V 


1 =5 

max 




0.25 



FIG. 7: The energy contained in the individual multipoles 
I — const (upper panel) as well as the energy radiated in all 
modes up to a given £ max (lower panel) as a function of the 
dimensionless mass parameter rj. 
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an artifact of the limited accuracy of the data and the 
polynomial fits. 

The significance of higher multipoles also reveals itself 
in the gravitational wave signal observed at fixed val- 
ues for the inclination angle 9 of the orbital plane of 
the binary. The combination of all multipoles in the 
signal ^4 shows a much more complex structure than 
the quadrupole on its own. We illustrate this in Fig. [5] 
where we plot real and imaginary part of 'J' 4 at r cx = 
36.5 Madm and 9 — 45.5°. The inclination angle has 
been chosen somewhat arbitrarily, but we get the same 
general result for other choices of 9. The pattern resulting 
from the inclusion of the octupole and higher order modes 
in the bottom panels of Fig. [8] shows a non-monotonic 
increase in the local maxima of the wave signal during 
inspiral which is remniscent of the pure quadrupole sig- 
nal in eccentric inspirals (cf. the "MayaKranc e02" entry 
in Fig. 1 of Ref. [lOdj |). In contrast there is little evidence 
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of eccentricity if we consider the individual modes of our 
simulation in Fig. [5J As a further impact of the higher 
order multipoles we note the significant increase in the 
amplitude between the top and bottom panel of Fig. [5J 
These observations qualitatively illustrate the impor- 
tance of including higher order multipoles in gravita- 
tional wave data analysis as has bee n po i nted out quan- 
titatively in the literature [H, S DM DM fiol Il04 
Il05l Il06l | 1 . In future work, we plan to investigate this is- 
sue more systematically using our numerical data in the 
context of matched filtering. 



IV. CONCLUSIONS 

In this work, we have pushed the mass ratio in numer- 
ical simulations of inspiralling black-hole binaries signifi- 
cantly beyond what has previously been published. Our 
simulations have been demonstrated to be convergent 
consistently with the discretization properties of the nu- 
merical code. We have resolved an issue specific to larger 
mass ratios by choosing a specific value of the shift damp- 
ing parameter. The overall errors in wave amplitude and 
phase as well as radiated energy and momenta are about 
5%. We have thus been able to validate existing fitting 
formulas for the amount of energy and linear momen- 
tum as well as the final spin of the merging hole in a 
range of the mass ratio previously unexplored. Within 
the error bars, we find our numerical results in agree- 
ment with phenomenological formulas. These are of high 
importance in the modeling of astrophysical phenomena, 
such as the growth of supermassive black holes. At least 
in the case of non-spinning holes, we conjecture that the 
fitting formulas can be used over the entire range of the 
dimcnsionless mass ratio parameter 77. 

Extending previous work [25| , we have further demon- 
strated that the percentage of gravitational wave energy 
radiated in higher order multipoles increases significantly 
as the mass ratio deviates from the equal-mass limit 
rj = 0.25 or q = 1. Quadratic fitting of our numeri- 
cal data indicates that about 32% of the total radiated 
energy will be contained in £ > 2 as the extreme mass 
ratio (EMR) limit is approached. For the q — 10 case at 
hand, we find about 25% of the energy to be contained in 
modes higher than the quadrupole. This distribution of 
the energy also manifests itself in the shape of the grav- 
itational waveform as measured for a fixed inclination of 
the binary orbit relative to the observer. The case 9 « 45 
degrees exhibits a significant change in the wave signal 
as we include higher order multipoles. While we only 
display results for one value of 9, we find this behavior 
to be similar for arbitrary alternative inclinations. 



FIG. 8: The Newman-Penrose scalar ^4 obtained at an in- 
clination of the orbital plane of 9 = 45.5 degrees including 
multipoles up to and including ^ m ax = 2 , 3, 4 and 5 (from 
top to bottom). 



Higher order harmonics as discussed in some of these references 
are not strictly equivalent to higher order multipoles but often 
imply the inclusion of additional multipoles 
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It will be important to extend the current study to 
spinning binaries in the future. This will allow us to 
address various important questions in astrophysics and 
gravitational wave physics. For example, there remain 
uncertainties about the magnitude of the recoil effect for 
spinning binaries of unequal mass [TEI, H3, [97|]. Also, it 
will be important to compare fitting formulas for final 
spin and recoil with numerical results for initially spin- 
ning black holes. 

The use of numerical waveforms in gravitational wave 
data analysis further requires careful comparisons with 
post-Newtonian predictions and the combinati on of nu- 
merical w ith p ost-Newtonian waveforms (e. g. [l05l Il07l 
HM EM fTTnt [TT1 HH Hi! . flT^ h In future work, we 
plan to perform such comparisons including various post- 
Newtonian techniques. We also believe that the current 
simulations facilitate approximate comparisons with per- 
turbative calculations of EMR binaries. Finally, an in- 
triguing question concerns the properties of black hole 
collisions or scattering experiments involving unequal- 
mass binaries traveling close to the speed of light. Stud- 
ies have so far been r estricted to collisions of nonspinning 
equal-mass binaries jl!5l . Ill6l ] and predict that as much 
as ~ 14% of the total energy of the system can be radi- 
ated in head-on collisions and even larger quantities for 
non-zero impact parameter. 

In summary, accurate simulations of black-hole bina- 
ries with mass ratio q = 10 are not only possible using 



current numerical techniques (if somewhat expensive), 
but also reveal a richness in structure beyond what has 
been observed in the nearly equal mass case. We con- 
sider our simulations of the non-spinning case to be the 
starting point of more exhaustive studies involving spins 
and/or eccentricity, which will be of significant value for 
current efforts to observe gravitational waves and im- 
prove our understanding of astrophysical questions in- 
volving the merger of black-hole binaries. 
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